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O : ABSTRACT 

Q ' Using high redshift radio sources as background, the 21 cm forest observations probe 

■ the neutral hydrogen absorption signatures of early structures along the lines of sight. 
Directly sensitive to the spin temperature of the hydrogen atoms, it complements the 
21 cm tomography observations, and provides information on the temperature as well 
as the ionization state of the intergalactic medium (IGM). We use a radiative transfer 
simulation to investigate the 21 cm forest signals during the epoch of reionization. We 
first confirmed that the optical depth and equivalent width (EW) are good represen- 

■ tations of the ionization and thermal state of the IGM. The features selected by their 

I ' relative optical depth are excellent tracers of ionization fields, and the features selected 

2 ■ by their absolute optical depth are very sensitive to the IGM temperature, so the IGM 

-4— > , 

c/3 , temperature information could potentially be extracted from 21 cm forest observation, 

■ 

thus breaking a degeneracy in 21 cm tomographic observation. With the EW statistics, 
^ we predict some observational consequences for 21 cm forest. From the distributions 

' of EWs and the number evolution of absorbers and leakers with different EWs, we see 

in 



O 

U 



o 



clearly the cosmological evolution of ionization state of the IGM. The number density 



CN ■ of potentially observable features decreases rapidly with increasing gas temperature. 

■ The sensitivity of the proposed EW statistic to the IGM temperature makes it a unique 



and potentially powerful probe of reionization. Missing small-scale structures, such as 



0^ , small filaments and minihalos that are unresolved in our current simulation, and lack 

' of an accurate calculation of the IGM temperature, however, likely have rendered the 

presented signals quantitatively inaccurate. Finally we discuss the requirements of the 
^ ' background radio sources for such observations, and find that signals with equivalent 

• widths larger than IkHz are hopeful to be detected. 
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1. Introduction 

The reionization of the intergalactic medium (I GM) is one of the landmark events in the history 



of the universe. While a lot of theoretical works (IWvithe &: Loeb 



Zhang et al.ll2007l:IChoudhurv fc Ferrarall2007l:IChoudhurv 



(|Fan et alJl2006l : lOunkley et alJl2009l : iGallerani et al 



2003 



Furlanetto et al 



2004a 



20091 . e.g.) and observational constraints 



200a, e.g.) h ave schemed up a globa l picture 



of the reionization process, many detailed issues remain uncertain (jCiardi &: Ferrarall2005l ). One of 
the most promising pr obes of the cosmic reionization is the redshifted 21 cm transition of HI (see e.g. 
Furlanetto et al.l 120061 for a review). The emission or absorption of the 21 cm radiation is directly 
related to the neutral component of the IGM at different redshifts, thus affording us the most clear 
view of the reionization history. Moreover, unlike the Lya resonance line which also traces the 
neutral hydrogen, the 21 cm line has an extremely small Einstein coefficient (^lo = 2.85 x lO^^^s"^), 
hence it is not saturated even at high redshift, and one can therefore safely predict the neutral 
fraction from its 21 cm optical depth. 

The most commonly discussed way of 21 cm observation is 21 cm tomography, in which the 
three dimensional structure of the en i ission or absorption of 21 era photons by the IGM against the 
CMB is observed dMadau et al.lll997l : ItozzI et ahlboool : llliev et ahlbood : iFurlanetto et al.ll2004al lb[). 



A less frequently discussed way of 21 cm observation is to observe the spectra of high redshift radio 
point sources, such as quasars or gamma ray burst (GRB) radio afterglows, and look for absorption 
or leak features produced by structures along the line of sight. Such 21 cm forest could be produced 
by structures of different scales, including the large scale HI/HII regions ("bubbles"), the cosmic 



the IGM at redshift z is (Field 



(Field 


1959: 


Madau et al. 


1997: 



Furlanetto et al.ll200fil ): 



hpC^Aio 



XHinuiz) 



327r kBvl Ts{l + z){dv\\/dr\\ 



(1) 



0.009(1 +5)(l + z)3/2^ 



H{z)/{l+z) 
dv\\/dr\\ 



where xhi is the neutral fraction of hydrogen at a specific position, nniz) and Ts are the number 
density of hydrogen and the spin temperature of the IGM at that position, respectively, and dun /dry 
is the gradient of the proper velocity along the line of sight, including both the Hubble expansion 
and the peculiar velocity. As we are going to study the large scale HI/HII regions, the peculiar 
velocity is neglected. In the second line, the factor (1 + 5) is the overdensity of the baryons, and 
Ts is in units of Kelvin. 

In the following, we will assume that the spin temperature always follows t he gas temperature , 



i.e. Ts = Tk- This is a good approximation f or the redshifts relevant here (ISantos et al 



where the Lya background from star formations (|Madau et al.lll997l : IChen &: Miralda-Escudell2004) ) 



200. 



is present, and the HI s pin temperature is s trongly cou pled to the kinetic temperature via the 
Wouthuysen-Field effect (jWouthuvsen Ill952l : lFieldlll959l ). The 21 cm forest signal is sensitive to 
the spin temperature, this is a great advantage compared with 21 cm tomographic observation 
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against the CMB background. The latter could also be sensitive to the spin temperature in the 
case of absorption, but once the spin temperature raises above the CMB temperature, the 21 cm 
signal become emission, which saturates and become almost independent of T5. This uncertainty 
in spin temperature ca uses large degeneracies between cosmological and astrophysical parameters 
in 21 cm tomography (jSantos &: Cooravl l2006l ) . The information extracted from the 21 cm forest 
observation may help us to determine the spin temperature at the relevant redshifts, and reduce 
the degeneracies in 21 cm tomography. 



Carilli et al.l (120021') used numerical sim ulation to illustrate the 21 cm forest imprinted on 



quasar spectrum. iFurlanetto &: Loebl (|2002l ) considered small scale features produced by mini- 



halos and galaxies. The absor ption features produced b y the HII regions ("ionized bubbles" 



before reionization was studied dFurlanettol bood : H bood ) using the analytical "bubble model" of 



reionization (^Furlanetto et al 



2004a 



in which the ionization proceeds inside-out, i.e. from high 



density regions to low density regions. 



The reionization process is very complicated, a.lthough the "bubble model" picture of reioniza- 
tion i s in general agreenient with simulatiori resul ts ( Zahn et al.ll2007 : Iliev et ahlboOT : Mellema et al 
20061 : iTrac fc CenI 120071 : IShin. Trac Cenll2008l ). the comphcated geometry of the ionized regions 
is not considered in this simple analytical model. Therefore, to extract information from future 21 
cm forest observations, it is imperative to examine how the 21 cm forest observables are correlated 
with the reionization process by using numerical simulations. In this paper, we study the large 
scale 21 cm forest produ ced by the neutral and ionized bubbles using the reionization simulation by 
Shin. Trac &: CenI ([200a). The high resolution simulation has detailed modeling of star formation 
and radiative transfer of ionizing photons. First, we check whether the equivalent width (EW) can 
be a representative indicator of the ionization and thermal state of the IGM, and to what extent the 
21 cm forests trace cosmic reionization. Then we use simulation data to predict the observational 
consequences for 21 cm forest. Comparing them with realistic parameters of radio arrays, we study 
what constraints we could put on reionization from 21 cm forest observations. 

We briefly describe the reionization simulation in section 2. In section 3, we check whether 
or not, and to what extent, the optical depth (hence the equivalent width) follows the ionization 
state Xi, density field p, and thermal history Tiqm of the IGM. From these analyses, we would see 
what information can be extracted from EW statistics. Then we predict the 21 cm forest signals - 
distributions and evolution of EWs - from the reionization simulation in section 4. In section 5 we 
investigate whether there are sufficient number of background sources that are luminous enough 
for such 21 cm forest observations. We summarize and discuss our results in section 6. 



Throughout the paper we use the sarne AC DM cosmology as IShin. Trac &: CenI (j2008l ) based 
on the WMAP3 results (see lSpergel et al.ll2007l . and references therein): {Qrn,^A,^b,h,as,ns) = 
(0.26, 0.74, 0.044, 0.72, 0.77, 0.95). 
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The Simulation 



In this paper, we use the reionization simulation described in detail in lShin. Trac &: CenI (|2008l ). 
The high-resolution simulation was run using a hybrid code that contains a N-body algorithm for 
dark matter, prescriptions for baryons and star formation, and a radiative transfer algorithm for 
ionizing photons. All components are solved simultaneously in a periodic box with comoving side 
length of 100 Mpc/h. 

The N-body calculations included 2880^ collisionless dark matter particles, each with a mass of 
3.02 X 10^ M^/h. This provided sufficient resolution to resolve halos down to masses of ~ 10*^ M^/h 
and account for the majority of photo-ionizing sources. For each particle, the local matter density 
and velocity dispersion are calculated in order to model the baryons and star formation. The 
baryons are assumed to have the same overdensity field as the dark matter, and the temperature 
field is set by the velocity dispersion of the particles. Radiative heating and cooling of the gas are 
not included in the hybrid simulation. 

The criterion for star formation was that the particles must have densities pm > WOpcritiz) 
and temperatures T > 10^ K. This cut in the temperature-density phase space effectively restricts 
star formation to regions within the virial radius of halos that cool efficiently through atomic line 
transitions. The simulation distinguishes between the first generation, Population III (Pop III) 
stars and the second generation. Population II (Pop II) stars. The stellar initial mass function 
is determined by following the metallicity enrichment of the halo and the IGM as described in 



Trac &: CenI (j2007l ). The input spectrum for ionizing radiation is divided in three energy ranges 
(eV): 13.61 < u < 24.59, 24.59 < u < 54.42, and u > 54.42, according to the ionization limits for 
for HI, Hel, and Hell. 

The radiative transfer algorithm for ionizing radiation is based on a phot on-advection schem e 



that is computationally much less expensive than the ray-tracing scheme in iTrac &: CenI (|2007l ) 



These two approac hes have been demonstra ted to give very similar results for the majority of the 



reionization epoch (jShin. Trac Sz Cenll2008l ). In the hybrid simulation, reionization is completed 
by z ~ 6. For post-analysis, the dark matter, baryons, and radiation are collected on a grid with 
720^ cells and the data is saved every 10 million years from z = 25 down to z = 5. In this paper, 
we make use of the gas density, ionization, and temperature fields. 

We mark a cell as "ionized" if the ionization fraction is greater than 0.5, and "neutral" other- 
wise. Slices of ionization field for several redshifts are shown in Fig. [H Each slice in the figure has 
a thickness of 139 kpc/h, and this applies to all the other slice-plots in the following. We see that 
the HII bubbles were highly biased and non-spherical; the universe became 50% ionized at redshift 
z ~ 9.5, and the non-gaussianity in the HII field was significant throughout the reionization. In ad- 
dition, the simulation also shows that Pop III galaxies played a very important role in the process, 
and a strong correlation existed between halo number density and bubble size for large bubbles. 



To characterize the ionized region statistically, we use the Priend-of-Friend (FoF) algorithm 
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z = 8.0 





z = 63 




Fig. 1. — Distributions of ionization fraction at redshifts z ~ 13.5, 11.2, 9.3, 8.0, 7.0, and 6.3. 
Highly ionized regions are represented by white while regions with ionization fraction below 50 
percent is shown as black. These are plots of slices in the simulation box with each side of 100 
Mpc/h. The global ionization fractions are ~ 0.094, 0.24, 0.46, 0.73, 0.91, and 0.98 respectively. 



to select ionized bubbles (introduced by llliev et al.l (j2006l )). which are made of connected ionized 
cells. The linking length is set to be 1.2 times the simulation data output cell size, so that cells 
are connected by one face at least. In addition, one bubble consists of 2 or more connected cells. 
Fig. [2] shows the size distributions of bubble s selected by ionized frac tion (xj-bubbles), with Xi > 
0.5, at several redshifts. As discovered by IShin. Trac &: CenI (j2008l ). before the percolation of 
the ionization bubbles, these ionized- fraction-selected bubbles (xj-bubbles) have the distributions 
with three distinct peaks. The characteristic volumes for the three peaks are 0.6, 0.03, and 0.006 
Mpc^/h^ respectively. These sizes do not change significantly with decreasing redshift, before one 
connected network of bubbles dominated at z < 10. The smallest peak in each distribution is due 
to limited resolution of the simulation (i.e. cell size) . The other two characteristic sizes are believed 
to be real features of the size distribution (jShin. Trac &: Cenll2008l ). 
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Fig. 2. — Number distributions of Xj-bubble volumes at redshifts z = 13.5, 11.2, 9.3, 8.0, 7.0, and 
6.3, respectively. 



3. The Equivalent Width as a Probe of Reionization 

The uniformly distributed IGM would produce only uniform absorption on the radio spectrum, 
it is the variation of density, ionization state, and spin temperature which produces varied absorp- 
tion that is observable. This is illustrated in Fig. O One may define an absorber or a leaker as a 
continuous part of the spectrum for which 

- rol > Tth (2) 

where tq is the average optical depth over the whole simulation box at that redshift, Tth is the 
threshold value. The strength of 21 cm absorption is typically not very strong due to its weak 
transition coefficient. We take r^h = 0.002, which could induce features observable with high signal 
to noise observation. The criteria — tq > Tth picks up absolute absorbers while the criteria 
To — Tu > Tth picks up absolute leakers. Alternatively, one can also use a relative definition 

'"^-"°'>i? (3) 



To 

where R is the threshold, and is taken as 0.5 here. We will see later (in Fig. [7]) that the morphology 
and distribution of the optically thin or thick regions remain nearly the same if we use a different 
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Fig. 3. — The optical depth as a function of observed frequency along a line of sight passing through 
the simulation box. The redshift here is 10.2. 



threshold. In this way, one can select relative absorbers and leakers. 

To characterize the strength of signals, we consider the equivalent width (EW) of both ab- 
sorbers or leakers. For absolute or relative absorbers, 

J absorber JIGM 

= I {l-e''°-^'')du (4) 

J absorber 

where ficAi is the residual flux of a high redshift radio source after the mean absorption by the 
IGM, and fi, is the flux at frequency v. For absolute or relative leakers, 

Wu = f {e^°~^^ - 1) dv. (5) 

J leaker 

Clearly ionized regions will produce leakers and neutral clumps will produce absorbers on high 
redshift quasar spectra. 

Before using the equivalent width as our statistics, we test whether it is a good tracer of the 
ionized (neutral) regions, and to what extent we can extract Xj, p, and Tjgm information from EW 





z = 63 




Fig. 4. — Distributions of relative optical depth at redshifts z ~ 13.5, 11.2, 9.3, 8.0, 7.0, and 
6.3. Regions with (tq — Ty)/TQ > 0.5 (optically thin) are represented by white while other regions 
are shown as black. These are plots of slices in the simulation box with each side of 100 Mpc/h. 
The global ionization fractions are same as Fig. [TJ The corresponding average optical depth for 
these redshifts are tq = 1.1 x 10"^ 7.5 x 10"^, 4.4 x 10"^, 1.8 x 10"^, 4.8 x 10"^, and 6.7 x 10"^ 
respectively. 



statistics. This test includes two parts: (1) the consistency between ionized regions and optically 
thin regions, and (2) the correspondence between 3D bubbles and the projected ID absorbers. 

To do the first test, we display slices of optical depth field in Fig. HI In the optical depth 
calculation, the spin temperature Tg is taken to be the temperature output from the simulation 
data. We have used the relative definition (tq — t)/to > 0.5 to define the optically thin cells. 
Optically thin regions are shown as white while optically thick regions are shown as black. We see 
good agreement between optically thin regions in Fig. H] and those highly ionized regions in Fig. [H 
and the configurations match very well between the Zj-slices and r-slices. For more quantitative 
comparison, we also made scatter plots of the absolute optical depth r and the relative optical 
depth (tq — t)/tq versus ionized fraction xi on a cell-by-cell basis in FigE] (for z = 9.3). There 
is a strong correlation between the relative optical depth and the ionized fraction, as well as a 
correlation between the absolute optical depth and the ionization fraction at Xj > 0.6. Note that 
cells with high optical depth (large lg(r) or very negative (tq — t)/tq) have low ionization fraction 
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Fig. 5. — Scatter plots of 21 cm optical depth r versus ionized fraction Xi (left panel) and the 
relative optical depth (tq — t)/tq versus ionized fraction x, (right panel). The redshift shown here 
is z = 9.3. 



Xi, while high threshold on (tq — t)/tq picks out mostly highly ionized regions. 

We may also use the FoF algorithm to select and connect cells by the value of the optical 
depth instead of the ionization fraction. Size distributions of the optically thin bubbles (denoted as 
relative r-bubbles) are shown in Fig. [S] for several redshifts. The size distribution of these relative 
T-bubbles are apparently very similar to that of the Xj-bubbles, but the peaks are more smoothed, 
probably due to the mixture of ionization state and density information in optical depth. So the 
optical depth in the relative definition is a reasonable representative of the ionization status of the 
IGM. Although it smears out the characteristic scales of ionized bubbles, the r-field does follow the 
large scale structure of the ionization field. 

We also tried different relative thresholds in selecting these r-bubbles, in order to see how the 
signals would change with the thresholds we used. Fig. \7\ shows the optical depth slices at redshift 
9.3 for thresholds (tq — t)/tq > 0.1, 0.3, 0.7, and 0.9, respectively. The shapes and distributions of 
the bubbles are quite similar among these slices with different relative optical depth thresholds. The 
T-slices with more strict thresholds have more dark dots (optically thicker regions) inside bubbles, 
and have less white dots (optically thinner regions) outside bubbles. All these r-bubbles match well 
with Xj-bubbles, and similar results can be found for other redshifts. So the observational features 
in 21 cm forests would not change significantly if different threshold is adopted. 

Features picked up with absolute threshold are more relevant to real observations. HIT Bubbles 
and HI islands selected by the absolute threshold \ti, — ro| > 0.002 for the corresponding 6 redshifts 
are shown in Fig. [8l Using this absolute optical depth threshold, there are both bubbles and islands 
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Fig. 6. — Number distributions of relative r-bubble volumes at redshifts z = 13.5, 11.2, 9.3, 8.0, 
7.0, and 6.3, respectively. The criteria to select these bubbles is (tq — t^)/tq > 0.5. 



at three higher redshifts (i.e. 13.5, 11.2, and 9.3), which generate leakers and absorbers on spectra 
respectively. At lower redshifts there are only islands, since the average optical depth itself is 
smaller than 0.002 (i.e. at redshifts 8.0, 7.0, and 6.3 here), so we could only observe absorbers. For 
the leakers, smaller threshold must be used in order to find them at these redshifts. The shapes 
and distributions of these r-bubbles are similar to those selected by relative thresholds. Combined 
with Fig. [1] and Fig. [H we come to the conclusion that no matter relative or absolute threshold 
we use, the optical thin regions follow the highly ionized region very well, at least on large scales. 
Size distributions of absolute r-bubbles and r-islands for the three higher redshifts and those of 
absolute r-islands for the other three lower redshifts are shown in Fig. [H 

In order to see whether the optical depth (hence the EW statistics) is a good tracer of density 
fluctuations, we plot density slices in Fig.llOi But because of the limited resolution of the simulation 
(the overdensity is smoothed within each cell), the density contrast is very small (typically 6 < 1) 
on these scales, and do not show much evolution, and there is little correlation with the ionization 
fraction. This means that we would miss strong absorbers produced by small but highly dense 
regions such as cosmic webs, this would affect the accuracy of this simulation in representing the 
real case. 
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Fig. 7. — Distributions of optical depth at redshift z = 9.3, when the average optical depth 
To = 4.4 X 10^^. Regions with relative optical depth (tq — t^)/to > threshold (optically thin) are 
represented by white while other regions are shown as black. The relative thresholds are 0.1, 0.3, 
0.7, and 0.9 from left to right and top to bottom, respectively. These are plots of slices in the 
simulation box with each side of 100 Mpc/h. 

Can we get some information about the temperature of the IGM from 21 cm forest signals? 
Although there is no heating information incorporated in the simulation, and the gas temperature is 
derived from the velocity dispersion of the particles, it is still possible to test the sensitivity of optical 
depth and the forest signal to the IGM temperature, by assuming several different gas temperatures. 
In top panels of Fig. \TT\ we plot relative r-slices at redshift 9.3 for the IGM temperatures of lO^K, 
the gas temperature from simulation, and lO^K, from left to right, respectively (in this excise we 
set the same temperature at every point within the simulation volume). Absolute r-slices with the 
corresponding gas temperatures are plotted in bottom panels of Fig. [TTl 

Obviously, r-slices with relative thresholds are insensitive to the gas temperature, while the 
absolute r-slices are very sensitive. When we raise the IGM temperature to be lO^K, the average 
optical depth at redshift 9.3 falls below the threshold 0.002, and no optical depth could drop below 
the mean value by 0.002, so there is no optically thin region in the absolute definition in this 
redshift, just as the three slices of lower redshifts in Fig. [HI If we raise the IGM temperature 
further to be lO^K, then very few regions could have optical depth larger than the average depth 
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z=135 z=11.2 z = 9.3 




Fig. 8. — Distributions of absolute optical depth at six redshifts: z = 13.5, z = 11.2, z = 9.3, 
z = 8.0, z = 7.0, and z = 6.3. Regions with tq — r^/ > 0.002 (optically thin) are represented by 
yellow, regions with Tu — tq > 0.002 (optically thick) are shown as dark blue, and the other regions 
lie in between are shown as light blue. The corresponding average optical depth for these redshifts 
are tq = 1.1 x 10"^, 7.5 x 10-^ 4.4 x lO'^, 1.8 x 10~^ 4.8 x 10"^, and 6.7 x 10"^ respectively. 
These are plots of slices in the simulation box with each side of 100 Mpc/h. 

by 0.002. Hence, we could easily extract Tk information from absolute r threshold observations 
while ionization state of the IGM could be extracted from relative r threshold observations. 

It is very interesting to find a method to measure the IGM temperature during the epoch 
of reionization. While the kinetic temperature of the IGM and spin temperature of hydrogen are 
usually assumed to be much higher than the CMB temperature in 21 cm power spectrum analysis, 
and the Tk information of the IGM is totally erased in "21 cm tomography" observations, the "21 
cm forest" observations serve as an excellent tool to separate the Tk information from the density 
and ionization status. 

Now we have seen some correspondence between the optical depth and the ionized fraction 
selected bubbles. However, in order to guarantee the validity of our EW statistics, we should also 
test the correspondence between observed absorbers and r-islands, or leakers and r-bubbles. To do 
so, we put random lines of sight through the simulation box, calculated the EW of each absolute 
absorber, and recorded their center positions. Then we checked whether the centers of absorbers 
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Fig. 9. — Size distributions of absolute r-islands at redshifts z = 13.5, 11.2, 9.3, 8.0, 7.0, and 6.3, 
respectively (blue solid lines), and those of absolute r-bubbles at the three higher redshifts z = 
13.5, 11.2, and 9.3 (red dot-dashed lines). The criteria to select these islands is — tq > 0.002, 
and the criteria for these bubbles is tq — Ti, > 0.002. 



reside in those r-islands we found. Results show that most EW features have their centers within 
r-islands, and those EW features which are not in r-islands are in one-cell- islands, which were 
not defined as islands in our "EOF" algorithm. Fig. 1121 shows scatter plots of absorbers' equivalent 
width and the volumes of their corresponding r-islands in which the centers of these absorb features 
reside. Absolute definition is used here for absorbers and islands, as is relevant in real observations. 
Although there is some correlation between the volume and equivalent width, we see that the scatter 
is very large, as the line of sight may pass through the region far from the center. One can also see 
a spike in the volume distribution, which corresponds to the connected HI region (partially ionized 
at lower redshifts) in the simulation box. This connected island does not disappear until the overall 
optical depth drops as seen at z = 7.0 when some places in the island failed to connect those regions 
around them and give birth to many small islands. More clearly, at z = 6.3, one connected island 
fragmented into three big islands and many smaller ones. They exist till late stages of reionization 
because of the self-shielding effect of these dense regions. 

In summary, optical depth can be a reasonable tracer of the ionized fraction and the equivalent 
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Fig. 10. — Distributions of overdensity at different redshifts. Overdense regions are represented by 
white while underdense regions are shown as black. 

width could serves as a great probe of cosmic reionization. While the features picked with relative 
thresholds arc better in extracting ionization information, the features picked by absolute thresholds 
are very sensitive to the IGM temperature. However, the small scale density fluctuation is not 
incorporated in our simulation, and the density field on large scales plays a minor role in determining 
the optical depth. So we do not expect to get information about local density fluctuation from EW 
statistics. In addition, we see correspondence between the 3-dimensional bubbles and the projected 
21 cm forest signals. 

4. Statistical Distribution of 21 cm Forest Signals 

First we calculate the average optical depth tq over the whole simulation box for each redshift 
z, and take this tq as the optical depth of the intergalactic medium, tjgm- Then we put 400 random 
lines of sights (LOS) through each side of the simulation box at a specific redshift z, and calculate 
the optical depths along each line of sight as a function of observed frequency. Using the threshold 
{tu — To) > 0.002, we pick every absorber along each line of sight, then calculate and record the 
EWs of them for every random LOS. The periodic boundary condition is employed here. 
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Fig. 11. — Distributions of optical depth with different gas temperatures at redshift z = 9.3. Top 
panels: optically thin regions with relative threshold {to — t,j)/tq > 0.5 are shown as yellow, optically 
thick regions with (ti, — tq)/tq > 0.5 are shown as dark blue, while other regions in between are 
labeled as light blue. The gas temperatures are assumed to be lO^K, gas temperature based on 
velocity dispersion, and lO^K, from left to right, respectively. Bottom panels: regions with absolute 
threshold tq — > 0.002 are shown as yellow, regions with — tq > 0.002 are shown as dark blue, 
while other regions are shown as light blue. The temperatures are the same as top panels. 



The number densities of absorbers and leakers are shown in Fig. [T3l The total number of 
these absorbers decreases significantly with decreasing redshift throughout the reionization. This 
is especially obvious for absorbers of small equivalent width, since the small islands were being 
ionized and becoming parts of ionized bubbles as the reionization proceeded. However, the number 
of absorbers with the largest equivalent width does not change as much, and the largest absorbers 
seem to persist to lowest redshifts. This is due to the self-shielding effect of the dense neutral regions, 
which precisely generated the spikes seen in Fig. [T2j They give us an opportunity to observe large 
absorbers even at lower redshifts. For the leakers, the number of smaller ones decreases with 
decreasing redshift since those small ionized bubbles were merging to generate larger ones, while 
the number density of larger leakers increases. 

Besides the equivalent width, the width of these absorbers or leakers is also an observable in 
radio probes, and it has good correspondence to the equivalent width as shown in Fig. [HI Just 
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Fig. 12. — Scatter plots of volumes of r-islands which hold absorbers versus EW of those absorbers 
at different reshifts. 



as EWs of absorbers and their corresponding r-islands were shown in Fig. [T2l Fig. [T3] shows the 
consistency between feature widths and the r-islands or r-bubbles associated with them. The 
maximum width of an absorber created by each island has a stronger correlation with the island 
size, and the same applies to relation between bubble size and the maximum width of leakers. 

The optical depth is very sensitive to the IGM temperature. We calculated the number density 
of absolute leakers that have Wiy > 1 kHz as a function of the gas temperature. The results for 
several redshifts are presented in Fig. [T5j The number densities of potentially observable absorbers 
decrease rapidly as the temperature increases. This means, on one hand, with the first generation 
of low frequency equipments, we may not be able to observe the signals that we are looking for 
if the IGM was really heated up, but on the other hand, we could obtain information on the gas 
temperature at the epoch of reionization from 21 cm forest observations. This is important because 
the IGM temperature contains very rich information about early structure formation, including the 
initial mass function of the first stars, properties of early mini-quasars and related X-ray heating, 
etc. In this sense, the 21 cm forest can be a powerful tool not only to probe the IGM in the early 
universe, but also to give constraints on the properties of early structures. 



Our prediction of the mean IGM optical depth is a little smaller than those by ICarilli et al. 

(j2002l ) at higher redshifts but larger than theirs at lower redshifts. For example, at z = 12 
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Fig. 13. — The number density of absolute absorbers (left panel) and leakers (right panel) with 
equivalent width larger than W^, as a function of W^. The curves for redshifts 13.5, 11.2, 9.3, 8.0, 
7.0, and 6.3 are represented by dotted (black) line, short-dashed (blue) line, solid (cyan) line, dot- 
dashed (red) line, long-dashed line (green), and long-and-short-dashed (purple) line respectively for 
absorbers, and only the three curves for higher redshifts are shown for leakers with the same line 
types as those for absorbers at the same redshifts. 



Carilli et al.l (|2002l ) found that tjgm ~ 1%, and some of the high-density regions reach r ~ 5% and 
above; by z ~ 8, the mean IGM optical depth has dropped to tjgm ~ 0.1%. In our simulation, 
however, tjgm ~ 0.9% at z ~ 12 and it only drops to t/ga/ ~ 0.18% at z = 8. 



However, unlike ICarilli et al.l (120021 ). we find very few strong absorbers. In our simulation 
box, only about one in two million cells w ould have r > 5% at z ~ 12, and no absorbers has 
r > 2% at z < 10, while ICarilli et al.l (j2002l ) still predicted 50 lines per unit redshift with r > 2% 
at z = 10 and about 4 lines per unit redshift at z = 8. Some of these differences may be due to 
the different values of cosmological parameters adopted in the two simul ation. Howev e r, the main 
difference probably comes from the different box size of the si niulations. ICarilli et al.l (|2002l ) used 
a simulation with smaller volume but higher spatial resolution (|Gnedinll200d ). in which the typical 
values of the cosmic overdensity 6 ~ 10. Our simulation has a much larger box volume, which is 
much more reliable for large scale reionization, but the size of each cell is larger. This affects the 
density fluctuation dramatically, and our typical overdensity within each cell is < 1, with a few 
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Fig. 14. — Top panels: the correspondence among the absorbers' width, their equivalent width, 
and the volumes of r-islands which created them. Bottom panels: the correspondence among the 
leakers' width, their equivalent width, and the volumes of r-bubbles which created them. The 
redshift shown here is 13.5. 
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Fig. 15. — Number density of absolute absorbers with equivalent width Wi, > 1 kHz as a function 
of the IGM temperature. The redshifts are z = 13.5 {dotted, black), z = 9.3 {dashed, blue), z = 
8.0 {.solid, cyan), and z = 7.0 {dot — dashed, red), respectively. 



regions have 6 > 1. As seen in Fig. \TU[ the density field evolves very little when smoothed within 
the cells, though at smaller scales there must be some evolution in density contrast. Therefore, the 
resultant local optical depth is underestimated. 

On the other hand, the first HII regions start to appear at reds hift about 8 in the simul ation 
by iGnedm hood), wh ich is much later than that in our simulation (IShin. Trac Cenll2008l ). As 



a result, ICarilli et al.l (120021) fo und higher optical depth of 21 cm forest than us at high redshift. 
Further, in ICarilli et al.l (|2002l ). those absorbers are typically filaments. They are at first shock- 
heated to about 100 K, and at z < 10 Lya heating increases the mean spin temperature to above 100 
K. In our simulation, the resolution is not high enough to capture the filaments anyway. The Lya 
heating is also unlikely to heat the gas to so high a temper ature, more likely it could only heat up the 



IGM to a temperature below 100 K before reionization (IChen fc Miralda-Escude 



2004 



Yueetal 



20091 ). In the simulation, these heating mechanisms are not incorporated, and the temperature is 
just from the velocity dispersion of the particles. Therefore we have lower spin temperature and 
this results in the higher mean IGM optical depth at lower redshifts, though we can artificially 
adjust the gas temperature as a proxy for incorporating the various heating mechanisms. 
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5. The number of high redshift radio sources 

To observe the absorption or leak features in the spectrum, the spectrum of the point source 
has to be observed with a certain precision, which depends on the brightness of the quasar. The 
minimum detectable flux density of an interferometer is related to the system temperature Tgyg, the 
effective aperture area ^effj channel width Afch, integration time t, and the signal-to-noise ratio 
S/N by: 

AC . - 2 kB Tsys ^ 

For the observation of resolved signals (i.e. the resolution pixel Az^ch is narrower than the width 
of the feature to be observed), the lower limit of the flux density of background sources for an 
observation of absorbers or leakers with |t;^ — tqI > Tf^ (when the t^, tq, and are all very small) 
is: 

1/2 



^cflf/^sys y V ^ 

where the ratio Agg/Tsys is an intrinsic parameter describing the sensitivity of an interferometry 
array. The systematic temperature is the sum of the receiver noise and the temperature of the 
uncleaned foregrounds. 

However, the optical depths of the "absorbers" and "leakers" are small, while the signals are 
relatively wide(10^ ~ 10^ kHz), so that observation of unresolved signals (or just below resolution) 
is more promising. In this case, the lower limit of the flux density for an observation of absorbers 
or leakers with > Wth is: 

1/2 



..„^n,w.(f^)(i^)(^) 

2 X lO^m^K-^X /lOOhry/^ 



(8) 

The demanded flux density of the background source is much lower compared to the resolved 
observations, making more high redshift radio sources qualified for such observations. 

We now investigate how many such quasars can be observed. Specifically, we calculate the 
number of quasars with flux density at the observed frequency = 1420.4/ (1 -|- z) MHz larger than 
the lower limit described above. In the following, we assume that all the radio quasars are of the 
same spectral energ y distribution with that of the powerful radio galaxy Cygnus A, i.e. oc u^^-^^ 



(jCarilli et al.ll2002l ). Then the limiting flux density converts to a limiting luminosity density at a 
rest-frame frequency of 151 MHz, P^i^- 



Jarvis et al.l (|200ll ) have pointed out that at z = 4, the comoving number density of radio 
sources with luminosity density at the frequency 151 MHz, Pi^i > 6 x 10^^ ergs s~^ Hz~^, is 2.4 x 
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10~^ Mpc~^. We extrapolate this number density to higher redshift and lower luminosity. Assuming 
that the comoving number density of qualified quasars at redshi ft z is given by nciy rf-Pi.^i > 
and according to the best-fitting luminosity function given bv lJarvis et al.l (|200ll ). we have 



-2.2 



(9) 



w here Cfz) is a norma lization constant, which in general is a function of redshift. But it is suggested 
bv lJarvis et al.l (|200ll ) that there is very little evidence for an abrupt decline in the comoving space 
density of the most luminous low-frequency-selected sources at high redshift, and the best-fitting 
model gives a constant comoving space density. If, optimistically, we hypothesize quasars evolve 
according to the fiat model, i.e. the number density does not vary with redshift in the comoving 
coordinate, then the number of qualified quasars (-P151 > -Pisi") in the whole sky per redshift 



interval is: 



dz 



(10) 



where rcuiz) = cdz/H{z) is the comoving coordinate c orresponding to redsh ift z. On the 
other hand, based on optically selected quasars from SDSS, ICirasuolo et al.l (j2006l ) suggest that 
independent of the adopted radio spectral index, the drop in the space density of optically bright 
radio loud quasars between 2; ~ 2 and z ~ 4.4 is at most a factor 1.5 — 2. If we adopt the upmost 
factor 2, then 



C{z) = C(z = 4)exp[-||(z-4)]. 



(11) 



We plot the number of quasars which satisfy the above criteria for probing signals with W^, > 1 
kHz in Fig.[T6]for the fiat evolution case in the left panel, and for the evolved case in the right panel. 
Here we choose three values of sensitivity, A^s/Tsys = 5000 m^K^^, 2000 m^K~^, and 1000 m^K~^ 
respectively, which can be achieved in future observation programs. The Square Kilometer Array 
(SKA) will have A^s/Tsys = 5000 m^K"^ at the frequency between 70 MHz and 300 MHz (see 
ihttp: / /www. skatelesc ope.org/} , and the Low Frequency Array (LOFAR) could marginally achieve 
Aes/Tsys = 1000 m^K"^ at 200 MHz (|http: / /www.lofar .org/j) . We can see that for the fiatly evolved 
quasar density, there would be many bright quasars which could be used as our background sources. 
However, for the steeply evolved quasar model, those suitable background quasars would be very 
few. 



6. Conclusions 



We used a numerical simulation with star formation and radiative transfer ( Shin. Trac &: Cen 



20081 ) to study the 21 cm forest signals during the epoch of reionization. We first defined absolute 
and relative absorbers and leakers on spectra of background sources produced by structures along 
the line of sight, and checked whether the absorbers or leakers can represent the ionization and 
thermal state of the IGM. We found that the optical depth can be a reasonable tracer of the ionized 
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Fig. 16. — The number of quasars in the whole sky that can be used to detect signals with 
VFjy > IkHz per redshift interval. Left panel: the number of qualified quasars in the flat evolution 
model. Right panel: the number in the steep evolution model. The sensitivity of the radio array is 
taken to be A^s/Tsys = 5000 m^K'^ 2000 m^K'^ and 1000 m^K~^ for the red-solid, blue-dashed, 
and black dotted line, respectively. The integration time is assumed to be 100 hours here. 



fraction of the IGM, and the equivalent width could serves as a great probe of cosmic reionization. 
The features selected by absolute threshold in optical depth are also very sensitive to the IGM 
temperature, so the temperature information of the IGM could be separated from the ionization 
fraction and density fluctuations with 21 cm forest observation, making it a good complement to 
21 cm tomography observation. At the scales we studied, the density fluctuation plays a minor role 
in determining the optical depth. Also, we see correspondence between the 3-dimensional bubbles 
and the projected 21 cm forest signals. 

The number densities of leakers and absorbers with different equivalent widths evolve with 
redshift, showing the evolution of ionization status of the IGM. While the total number of absorbers 
decreases significantly with decreasing redshift, the largest absorbers persist to low redshifts because 
of the self-shielding effect. As for the leakers, diminished small ones and ballooning large ones with 
time show the signature of merging ionized bubbles. From the largest width of leakers or absorbers 
in 21 cm forest spectra, we may infer the size of ionized bubbles or neutral islands at specific 
redshifts. 
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The most important advantage regarding the 21 cm forest observation is its sensitive depen- 
dence of signals on the IGM temperature. The number density of potentiahy observable signals 
decreases dramatically with the increasing gas temperature, so the temperature of the IGM at each 
redshift could be constrained potentially through the 21 cm forest observations. As the thermal 
history of the IGM carries rich information about the first luminous sources, including the popula- 
tion of the first stars, properties of mini-quasars, and complex heating processes, the 21 cm forest 
observation is a potentially powerful probe of the early structure formation of the universe. The 
measurement of the gas temperature would also help improve the precision of the measurement of 
cosmological parameters. 



We also compared our results with a previou s work by ICarilli et al.l (|2002l ). who used a simu- 
lation with higher resolution but smaller volume (jGnedinll200Cl ). While our simulation has a much 
larger volume which is necessary to account for large HII regions at late epochs of reionization, 
small-scale structures such as small filaments and minihalos are not resolved, especially in the 
binned data used for post-analysis. With higher resolution simulations, higher gas overdensities 
and hence higher optical depths would be expected. Since heating processes are not included in 
our simulation, the assumed gas temperature is not realistic, but we do consider a range of possible 
temperatures. With these limitations, we can not make an accurate prediction for the signal, but 
instead provide a basic expected range. Nonetheless, our simulation does reveal some qualitative 
features of the large scale 21cm signal, and propose here that the 21 cm forest observations, es- 
pecially the EW statistics, could play an important role in extracting temperature information of 
the IGM. We reserve more realistic investigation of the 21 cm forest to future works, for exam- 
ple by using higher-reso lution versions of the hydro-|-radiative transfer simulations presented in 



Trac. Cen &: Loebl (120081 ). which have self-consistent hydrodynamics and temperature evolution. 



Finally we discussed the the number of luminous high redshift radio sources which is potentially 
available for observing such forest signals. We investigated the requirement on the luminosity of 
background radio sources for the observation of signals with equivalent widths larger than IkHz. 
We considered how the number of qualified quasars in the whole sky evolves with redshift according 
to two models. If quasars evolves according to the flat evolution model, we would have sufficient 
quasars as background sources at high redshift; if quasars evolve fast, however, only those signals 
at relatively lower redshift can be observed in the near future. 



We thank Bin Yue, Andrea Ferrara, Steven Furlanetto and Chris Carilli for helpful discussion. 
This research is supported in part by the NSFC grant 10525314, 10503010, 10373001, 10773001, 
by the CAS grant KJCX3-SYW-N2, MoST 973 grant 2007CB815401, and by NASA ATP grant 
NNG06GI09G. Computing resources were in part provided by the NASA High-End Computing 
(HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research 
Center. 



- 24 - 



REFERENCES 

Carilli, C. L., Gnedin, N. Y., & Owen, F. 2002, ApJ, 577, 22 
Chen, X. & Miralda-Escude, J., 2004, ApJ 602, 1 
Choudhury, T. R. 2009, larXiv:0904.4596l 
Choudhury, T. R. & Ferrara, A. 2007, MNRAS, 380, L6 

Ciardi, B. & Ferrara, A. 2005, SSRv, 116, 625 (updated in 2008: astrp-ph/040901 8'^^2) 

Cirasuolo, M., Magliocchetti, M., Gentile, G., Celotti, A., Cristiani, S., & Danese, L. 2006, MNRAS, 
371, 695 

Dunkley J. et al. 2009, ApJS, 180, 306 

Fan, X., Carilli, C, & Keating, B. 2006, ARAA, 44, 415 

Field, G. B. 1959, ApJ, 129, 525 

Furlanetto, S., 2006, MNRAS, 370, 1867. 

Furlanetto, S. R. & Loeb, A. 2002, ApJ, 579, 1 

Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, PhR, 433, 181 

Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004a, ApJ, 613, 1 

Furlanetto, S. R., Zaldarriaga, M., &; Hernquist, L. 2004b, ApJ, 613, 16 

Gallerani, S., Salvaterra, R., Ferrara, A., and Choudhury, T. R. 2008, MNRAS, 388, L84 

Gnedin, N. Y. 2000, ApJ, 535, 530 

Iliev, I. T., Mellema, G., Pen, U.-L., Merz, H., Shapiro, P. R., k Alvarez, M. A. 2006, MNRAS, 
369, 1625 

Iliev, I. T., Mellema, G., Shapiro, P. R., & Pen, U.-L. 2007, MNRAS, 376, 534 

Iliev, I. T., Shapiro, P. R., Ferrara, A., & Martel, H. 2002, ApJ, 572, L123 ApJ, 619, 684 

Jarvis, M. J., Rawlings, S., Willott, C. J., Blundell, K. M., Eales, S., k Lacy, M. 2001, MNRAS, 
327, 907 Dutta, S. 2007, ApJ, 670, 39 

Madau, P., Meiksin, A., k Rees, M. J. 1997, ApJ, 475, 429 

Mellema G., Iliev I. T., Pen U. L., k Shapiro P. R. 2006, MNRAS, 372, 679 

Mesinger, A. k Furlanetto, S. 2008a, MNRAS, 386, 1990 



- 25 - 



Miralda-Escude, J., Haehnelt, M., & Rees, M. J. 2000, ApJ, 530, 1 

Santos, M. G., Amblard, A., Pritchard, J., Trac, H., Cen, R., & Cooray, A. 2008, ApJ, 689, 1 
Santos, M. G. & Cooray, A., 2006, PRD 74, 083517. 
Shin, M.-S., Trac, H., & Cen, R. 2008, ApJ, 681, 756 

Spergel, D. N. et al. 2007, ApJS, 170, 377 Ohta, K., k Hattori T. 2006, PASJ, 58, 485 

Tozzi, P., Madau, P., Meiksin, A., & Rees, M. J. 2000, ApJ, 528, 597 

Trac, H. & Cen, R. 2007, ApJ, 671, 1 

Trac, H., Cen, R., & Loeb, A. 2008, ApJ, 689, L81 

Wouthuysen, S. A., 1952, AJ, 57, 31. 

Wyithe, J. S. B. &Loeb, A., 2003, ApJ, 586, 693. 

Xu, Y., 2006, "21cm forest" (in Chinese), B.Sc. thesis of Peking University, supervised by X. Chen 
and Z. Fan. 

Yue, B., Ciardi, B., Scannapieco, E., Chen, X., 2009, "21 cm signal from the IGM and minihalos", 
lar Xiv: 0906 . 3 1 05 [astro-ph. C O] , MNRAS accepted 

Zahn, O., Lidz, A., McQuinn, M., Dutta, S., Hernquist, L., Zaldarriaga, M., & Furlanetto, S. R. 
2007, ApJ, 654, 12 

Zhang, J., Hui, L., & Haiman, Z. 2007, MNRAS, 375, 324 Int. J. Mod. Phys. A, 19, 2385 



This preprint was prepared with the AAS lATp^X macros v5.2. 



